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Abstract 

The  Tnertial  Upper  Stage  (IUS)  being  developed  for  use 
aboard  the  Space  Shuttle  is  composed  of  three  solid  fuel 
stages  plus  a  satellite  payload.  One  mission  of  the  IUS 
system  is  to  launch  from  a  shuttle  parking  orbit  and  place 
the  satellite  in  geosynchronous  orbit  in  minimum  time.  Actual 
Space  Shuttle  parking  orbit  data  and  IUS  characteristics  were 
used  in  this  study  to  examine  the  sequential  timing  and  orien¬ 
tation  in  inertial  space  of  each  stage  as  it  is  fired  while 
the  spacecraft  moves  along  a  transfer  orbit  to  geosynchronous 
orbit.  In  addition,  the  sensitivity  of  the  total  transfer 
time  and  the  final  orbital  state  was  found  as  a  result  of  not 
meeting  one  or  all  of  the  time  and  orientation  parameters. 

This  problem  is  unique  in  that  it  considers  an  optimal 
orbit  transfer  problem  involving  solid  fuel  stages  of  fixed 
thrust  and  burn  time.  Previous  work  with  liquid  fuel  engines 
examined  orbital  transfers  with  the  intent  of  minimizing  the 
amount  of  propellant  or  required  velocity  change  needed  to 
accomplish  the  transfer. 


THREE  RURN  TNF.RTTAL  UPPER  STAGE 
OPTIMAL  ORBIT  TRANSFER 

I .  Introduction 

Previous  work  involving  optimal  orbit  transfers  was 
centered  around  liquid  fuel  engines.  The  driving  concern 
with  a  liquid  system  is  to  minimize  the  amount  of  propellant 
to  be  carried  on  the  spacecraft  to  keep  the  overall  weight  as 
low  as  possible.  This  is  especially  important  for  spacecraft 
of  one  stage,  where  the  engine  is  required  to  be  restarted 
for  additional  burns.  There  are  a  number  of  studies,  of 
which  Reference  9  is  an  example,  that  explore  the  transfer 
from  one  orbit  to  another  using  two  impulsive  burns.  One 
study  was  found  which  considered  three  burn  optimal  trans¬ 
fers  (Ref  10),  but  minimized  the  total  velocity  change  needed. 

The  advent  of  the  Space  Shuttle  made  it  possible  to 
carry  a  small  spacecraft  into  near  earth  orbit  inside  the 
cargo  bay  and  then  launch  it  from  that  parking  orbit  into 
another  orbit.  The  storage  problems  associated  with  liquid 
fuel  systems  made  solid  fuel  engines  preferable  for  this 
type  of  spacecraft.  The  initial  configuration  for  the 
Inertial  Upper  Stage  consisted  of  two  solid  fuel  engines. 

This  allowed  for  the  use  of  previous  orbital  transfer  studies 
to  define  the  exact  capabilities  of  the  configuration.  The 
desire  to  expand  this  capability  led  to  the  design  of  the 


1 


three  stage  IU5  which  necessitated  the  need  to  examine  a 
three  burn  optimal  transfer  to  geosynchronous  orbit  in  mini¬ 
mum  time. 

As  in  past  studies,  the  firing  of  each  stage  of  the 
spacecraft  is  treated  as  an  impulsive  burn.  This  assumption 
is  justified  because  the  burn  time  of  a  rocket  engine  is  very 
small  when  compared  to  the  total  time  spent  in  an  orbital 
transfer.  In  addition,  transfers  in  the  near  vicinity  of  the 
earth  are  normally  analyzed  in  three  phases  which  yield  in¬ 
creasing  degrees  of  accuracy. 

The  first  phase  is  to  consider  the  earth  and  the  space¬ 
craft  as  being  alone  in  inertial  space.  The  earth  is  con¬ 
sidered  perfectly  spherical  and  homogeneous,  and  the  effects 
of  the  other  major  bodies  such  as  the  sun  and  moon  are  ignored. 
These  assumptions  allow  for  a  two-body  analysis  of  an  object 
under  the  influence  of  a  single  gravitational  field.  The 
equations  of  motion  used  to  describe  the  behavior  of  an 
object  in  this  environment  are  accurate  enough  to  allow  for 
initial  mission  planning.  The  second  phase  is  to  take  the 
two-body  equations  of  motion  and  add  the  oblateness  and  non- 
homogeneity  effects  on  the  earth.  The  results  arc  more 
accurate,  but  they  differ  only  slightly  from  those  obtained 
from  the  two-body  equations  alone.  The  last  step  is  to  again 
modify  the  equations  of  motion  to  account  for  the  effects  of 
the  sun,  moon,  and  other  major  bodies.  The  end  product  i  , 
indeed,  accurate  but  extremely  expensive  and  tedious  to  extract. 

This  study  follows  the  tradition  of  previous  work  by 
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assuming  impulsive  burns  and  two-body  orbital  behavior  Tor 
this  initial  mission  analysis. 
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II.  Problem  Statement 


The  starting  place  for  this  problem  is  the  near  earth 
parking  orbit  of  the  Space  Shuttle.  The  IUS  System  Program 
office  at  Patrick  AFB  (Ref  1)  was  kind  enough  to  provide  the 
following  data  describing  the  shuttle  parking  orbit  at  the 
moment  of  orbit  attainment: 


TABLE  I 

Initial  Parking  Orbit  Data 


latitude 

21.3624° 

longitude 

59.2825° 

altitude 

7.11693  x  105  ft 

inclination 

28.78871° 

eccentricity 

0.01497 

velocity 

2.561851  x  104  ft/sec 

x  =  1.029312  x  107  ft 

xd  =  -2.248185  x  104  ft/scc 

y  =  1.732354  x  107  ft 

yd  =  9.356206  x  103  ft/sec 

z  =  7.881747  x  106  ft 

zd  =  7.958385  x  103  ft/sec 

where  x,  y,  z,  xd,  yd,  and  zd  are  the  position  and  velocity 
components  in  the  earth- centered  inertial  (ECI)  coordinate 
frame.  In  order  to  describe  subsequent  positions  and  times, 
this  initial  fix  (i.e.  initial  conditions)  was  chosen  as  a 
starting  place  for  the  orbital  problem  at  time  equal  to  zero. 
For  the  IUS,  the  engine  burn  times  are  of  fixed  duration  as 
shown  in  Table  IT;  therefore,  the  time  of  burn  initiation  and 
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TABU:  TT 


IUS  Characteristics 

Vehicle  Specific  Burn  Propellant 


Stage  Weight  Impulse  Time  Delta-V  Weight 

_ (lbs) _ (sec) _ (sec)  (ft/sec) _ ( lbs) 


1,2,  fi 3 

56119 

294.332 

123.91 

4242.175 

20263 

2$  3 

33813.6 

294.189 

140.35 

9565.712 

21505.6 

3 

8887. 5 

302.492 

94.49 

10997.98 

6016.6 

the  orientation  of  each  burn  in  inertial  space  are  the  vari¬ 
ables  of  interest. 

This  problem  is  worked  entirely  in  an  ECI  frame  of 
reference.  Figure  1  describes  the  orientation  of  each  stage 
prior  to  engine  initiation.  Since  the  assumption  has  been 
made  that  each  burn  will  he  treated  as  impulsive,  this  corro 
ponds  to  orienting  the  Delta-V  vector  in  inertial  space. 


z 


Fig  1.  Spherical  Coordinates  Definition 
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The  angles  alpha  and  beta  shown  in  Fig  1  completely  describe 
the  orientation  of  the  spacecraft  and  its  Delta-V  vectors  in 
inertial  space.  A  third  parameter,  time,  can  then  be  used  to 
pinpoint  the  position  of  the  spacecraft  along  each  segment  of 
the  transfer  orbit.  A  pictorial  representation  of  this  is  as 
follows : 


Fig  2.  Transfer  Orbit 


where 

T1  =  coast  time  to  firing  of  first  stage 
A1  5  B1  =  orientation  of  first  stage  at  firing 
T2  =  coast  time  to  firing  of  second  stage 
A2  f?  B2  =  orientation  of  second  stage  at  firing 
T3  =  coast  time  to  firing  of  third  stage 
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A3  f|  R3  =  orientation  of  third  stage  at  firing 
AV1,AV2,AV3  =  impulsive  velocity  changes  per  burn 

This  allows  the  entire  orbit  from  initial  parking  orbit  to 
final  position  to  be  described  by  the  nine  parameters  T1 ,  A1 , 
Bl,  T2 ,  A2,  B2 ,  T3 ,  A3,  and  B3. 

As  stated  before,  the  calculation  of  the  satellite's 
position  and  velocity  as  it  moves  through  each  segment  of  the 
transfer  orbit  will  be  done  using  a  two-body  formulation. 
Several  prediction  algorithms  are  in  existance  today;  all  of 
them  solve  variations  of  the  differential  equation 

R  *  -  -Hr  R  (1) 

K 

where  v  is  the  geocentric  gravitational  parameter,  TT  is  the 
vector  position  of  the  satellite,  R  is  the  vector  accelera- 
of  the  satellite,  and  R  is  the  scalar  distance  from  the  center 
of  the  attracting  body  to  the  satellite.  Two  such  algorithms 
were  used  in  this  study.  The  first  is  an  actual  series  solu¬ 
tion  of  the  two-body  differential  equations  known  as  the  f 
and  g  series  solution  to  the  Kepler  problem  (Ref  3").  The  f 
and  g  scries  solution  affords  a  rapid  means  of  determining, 
position  and  velocity  as  a  function  of  time  and  just  as 
importantly,  avoids  the  penalties  associated  with  the  numeri¬ 
cal  integration  of  differential  equations.  A  second  prediction 
algorithm,  known  as  Hamilton's  two  body  equations  (Ref  2),  is 
composed  of  six  first-order  differential  equations  and  was 
used  only  to  check  the  accuracy  of  the  subroutine  written  to 
do  the  prediction  problem. 
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By  using  the  f  and  g  scries  in  a  subroutine,  an  algor¬ 
ithm  defining  satellite  position  and  velocity  as  a  function 
of  all  nine  parameters  can  be  formed. 

Algorithm 

1.  Given:  initial  position,  velocity,  and  time  (X 
and  TO  =  0) 

2.  Coast  in  parking  orbit  for  time  T1  and  calculate 
final  position  and  velocity. 

3.  Align  the  first  burn  in  inertial  space  using  A1 
and  B1  and  calculate  its  component  velocity  con¬ 
tribution.  Add  this  to  the  velocity  from  (2)  to 
obtain  a  new  state  vector. 

4.  Coast  in  new  orbit  for  time  T2  and  calculate  a 
final  state  vector.  Repeat  the  procedure?  used 

in  (3)  for  the  second  stage  and  find  the  new  orbit 
parameters . 

5.  Coast  in  final  orbit  for  time  T3,  add  the  vectornl 
velocity  contribution  of  the  third  stage,  and  cal¬ 
culate  the  final  position  and  velocity. 

The  requirement  that  the  satellite's  final  orbit  be  geosyn¬ 
chronous  yields  five  constraints  on  the  final  state  vector 
as  defined  in  the  P.CI  frame. 

1.  The  z  component  of  position  must  equal  zero. 

2.  The  zd  component  of  velocity  must  equal  zero. 

3.  The  satellite's  speed  must  equal  geosynchronous 
speed . 

4.  The  satellite's  distance  must  equal  geosynchronous 
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distance. 

5.  The  final  orbit  must  be  circular. 

The  problem  then  is  to  find  values  for  the  parameters  T1 ,  Al, 
Bl,  T2,  A2,  B2,  T3,  A3,  and  B3  that  allow  the  satellite  to 
transfer  from  its  initial  position  in  the  parking  orbit  to  a 
final  orbit  which  satisfies  the  specified  end  conditions  in 
minimum  time.  It  is  readily  apparent  that  this  problem  is 
not  amenable  to  closed  form  analytic  treatment,  which  implies 
the  need  to  solve  it  iteratively. 


III.  The  Optimization  Problem 

A  good  definition  of  a  parameter  optimization  problem 
is  how  to  change  parameters  in  order  to  satisfy  end  conditions 
with  the  least  effort.  Optimization  theory  provides  the  dir¬ 
ection  needed  to  change  each  parameter  to  reduce  end  constraint 
errors  and  also  offers  figures  of  merit  that  can  be  used  to 
tell  how  well  the  procedure  is  progressing.  This  problem,  as 
stated,  involves  finding  the  values  of  nine  parameters  that 
force  the  final  orbit  to  satisfy  five  equality  constraints 
and  allow  the  satellite  to  complete  the  orbital  transfer  in 
minimum  time.  This  defines  the  performance  index  as 

G  =  T1  +  T2  +  T3  (2) 

which  is  the  total  time  of  flight.  The  initial  conditions 
for  the  problem  are  simply  the  initial  position  and  velocity 
components  supplied  for  the  parking  orbit  for  reference  time 
equal  to  zero.  The  final  conditions  or  end  constraints  are 
labeled  by  the  vector  M,  where 

M( 1)  =  zd 
M( 2)  =  z 

M(  3)  =  (xd2  +  yd2  +  zd2)  1/2-  1.0096  x  1 04  (  ^ 

M(4)  =  (x2  ♦  y2  +  z2)1/2-  1.3811  x  108 
M(5)  *  x(xd)  +  y(yd)  +  z(zd) 

g 

where  1.3811  x  10  feet  is  geosynchronous  orbital  distance. 
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and  1.0006  x  10^  feet  per  second  is  synchronous  speed.  F.ach 
constraint  is  written  so  that  when  its  numerical  value  equals 
zero,  the  desired  condition  has  been  obtained.  By  calculating 
the  norm  of  the  M  vector,  a  figure  of  merit  is  produced  which 
indicates  the  degree  to  which  the  end  constraints  arc  satis¬ 
fied.  The  optimization  problem  is  then  to  drive  the  norm  of 
the  M  vector  to  zero  and  the  performance  index  to  its  minimum 
value. 

One  procedure  for  accomplishing  this  is  referred  to  as 
suboptimal  control.  In  an  optimal  control  problem,  the  desired 
controls,  which  are  functions  of  time,  are  directly  calculated. 
For  this  problem  that  would  mean  the  alpha  and  beta  parameters 
are  functions  of  time.  The  assumption  that  all  three  burns 
are  impulsive  converts  the  alpha  and  beta  controls  to  scalar 
parameters.  The  suboptimal  control  technique  approximates 
these  controls  using  polynomials. 

In  this  problem  the  controls  are  the  three  pairs  of 
angles  needed  to  align  the  stages  for  firing  and  the  times 
of  firing.  If  all  the  parameters  are  expressed  in  one  vector 
A,  then 

A  =  [T1  A 1  B1  T2  A2  B2  T3  A3  B3]T  (4) 

This  vector  A  then  contains  all  the  information  needed  to 
find  the  final  position  and  velocity  of  the  satellite. 

Xf  =  X  f(  A)  (5) 

In  addition,  the  vector  A  also  supplies  all  the  information 
needed  to  evaluate  the  performance  index  and  the  end  condition 
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constraints . 


G  =  C(A) 
M  =  M(A) 


(6) 


If  an  augmented  performance  index,  F,  is  defined  functionally 
as 

F(A,v)  =  G  (A)  +  vTM(A)  (7) 

where  v  represents  a  vector  of  Lagrange  multipliers,  it  must 
satisfy  the  first  variational  requirements 


faT(A,v)  =  0 


..  3F 
A  9A 


FvT(A,v)  =  M(A)  =  0 

to  he  a  minimum  solution.  F^  is  the  partial  of  the  augmented 
performance  index  with  respect  to  the  parameter  vector  A. 

The  F^  matrix  represents  the  change  in  performance  resulting 
from  a  change  in  each  parameter  in  A  (gr  ad  i  ent ")  .  Since  there 
are  nine  parameters  in  A,  is  a  row  vector  of  nine  elements 
Hull  and  F.dgeman  (Ref  6")  describe  a  second-order  parameter 
optimization  technique  and  algorithm  specifically  for  appli¬ 
cation  to  suboptimal  control  problems.  Johnson  (Ref  7)  and 
Peterson  (Ref  8)  used  this  algorithm  to  develop  a  computer 
program  to  analyze  aircraft  time  to  turn  problems.  Their 
program  uses  three  controls,  differential  equations  of  motion 
and  two  end  condition  constraints.  The  program  was  modified 
to  accommodate  nine  parameter  controls,  a  scries  solution  to 
differential  equations  of  motion,  and  five  end  condition 


constraints.  In  short,  the  program  which  adopts  Hull  and 
Edgeman's  algorithm  uses  second-order  information  to  deter¬ 
mine  how  to  change  the  parameters  in  the  vector  A  and  the 
Lagrange  multiplier  vector,  v,  such  that  the  vectors  F^  and 
M  are  driven  to  zero.  The  vector  relationships  used  to 
change  A  and  v  are  derived  from  Eq  (8).  For  any  guessed 
values  of  A  and  v  that  are  not  a  solution,  the  equality  will 
not  hold. 

FAT  *  0  M(A)  t  0  (9) 

If  Eq  (9)  is  then  linearized  about  A  and  v,  then 


and  define 


iFAT  *  FAA5A  *  "*T 

6M  =  M.  6  A 
A 


«PAT  -  -  Pr/ 


<5M  =  -QM 

where  P  and  Q  are  scaling  factors,  yields 

FaaSA  *  mat«v  -  -pf/ 
MaSA  =  -QM 


which  can  he  solved  for  fiv  and  <SA. 

-  (MArAA-1MAT)‘1(-PMAPAA-,VAT  .  QM) 

*A  •  "faa"1(I'faT  4  V  ^ 


M 


3M 
A  TK 


and 
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AA 


92F 

32A 


(10) 

(ID 

(12) 

(13) 

(Id) 

(15) 


(16) 

(17) 


where 


and  P  and  Q  are  scaling  factors  which  control  optimization 
and  end  condition  satisfaction.  FAA  represents  a  change  in 
slope  and  the  direction  in  which  it  is  increasing  or  de¬ 
creasing.  The  algorithm  developed  to  iteratively  change  A 
and  v  follows: 

1.  Guess  A  and  v 

2.  Use  the  f  and  g  series  solution  to  find 

3.  Compute  M,  MA,  M^,  FA,  and  FAA 

4.  Pick  values  for  P  and  Q  and  compute  6v  and  6A 

5.  Set  A  =  A  +  6A  and  v  =  v  +  <Sv 

6.  Check  convergence  and  if  unsatisfied,  go  to 
step  (2) 

The  computer  program  uses  this  algorithm  plus  a  gradient 
approach  that  allows  the  direct  calculation  of  an  initial  v 
vector  which  eliminates  guessing  five  parameters.  The  pro¬ 
cedure  for  doing  this  was  also  developed  by  Hull  and  F.dgeman 
and  uses  first-order  information  as  follows: 


v  -  (MaMaT)'1[(Q/P)M  -  MaGaT]  (18) 


where 


„  _  3C. 

CA  9  A 


(19) 


For  a  typical  problem,  the  gradient  portion  of  the 
program  would  be  used  to  locate  the  vicinity  of  the  functional 
minimum.  At  this  point  where  gradient  methods  lose  their 
effectiveness,  the  second-order  algorithm  would  be  used  to 
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rapidly  converge  the  problem.  As  this  procedure  is  followed, 
the  scaling  factors  P  and  Q  are  gradually  raised  from  very 
small  values  to  final  values  of  one.  Tn  the  final  convergence 
cycles,  P  and  Q  must  be  equal  to  one  to  yield  an  optimal  solu¬ 
tion.  Although  this  procedure  would  appear  straightforward, 
its  application  to  specific  problems  can,  at  times,  require 
some  finesse  as  will  be  shown  later. 


Numerical  Methods 
'(Kef  6:484) - 


The  first-  and  second-order  optimization  routines  re¬ 
quire  the  matrices  M,  M^,  M^A,  and  F^.  These  matrices  were 
evaluated  using  numerical  techniques.  F^  and  are  defined 

as 


F 


A 


GA  + 


(20) 


AA 


AA 


5 

.  ViMAA- 

i  =  l  i 


(21) 


since 


and 


G  =  T1  +  T2  +  T3 

ga  =  [1  0  0  1  0  0  1  0  0] 

GAA  =  ^ 0  ^  9x9 


(22) 

(23) 

(24) 


The  only  quantities  that  need  to  be  calculated  are  M,  M^, 
and  MAA- 

The  M  matrix,  or  error  matrix,  is  easily  evaluated 
after  determining  the  final  position  and  velocity  vectors 
using  the  f  and  g  scries  solution  to  the  Kepler  problem. 


The  and  matrices  arc  determined  using  a  central 


differences  numerical  derivative  technique.  The  technique 


uses  the  initial  values  of  the  A  vector  parameters  (A  )  to 
calculate  an  initial  M.  P.ach  A^  is  then  positively  perturbed 


by 


A  .  =  A  +  6 
n+  n  n 


(25) 


and  then  negatively  by 


n- 


A  -  6 
n  n 


(26) 


Using  An+  and  A  ,  M+  and  M_  are  calculated.  The  central 
differences  representation  for  is  then  given  by 


M 


A 


n 


n 

M  -  M  , 

*  '  +  «(«„) 
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(27) 


n 


2  2 
where  o(6  )  represents  an  error  term  of  order  6^  where  6 

v  n  n  n 

is  a  very  small  positive  number.  The  matrix  contains  five 
rows.  The  first  row  is  determined  using  in  Eq  (27),  the 
second  row  using  and  so  on.  The  matrices  are  deter¬ 

mined  in  a  similar  manner;  however,  two  elements  A^  and  A^ , 
must  be  perturbed  both  positively  and  negatively  to  obtain 
M++,  M+  ,  M  ,  and  M  _.  The  central  differences  representa¬ 


tion  for  is 

n"  m 


M 


A  A 
n  m 


M  -  2M  +  M  7 

-  "  T-  ~  *  > 

n 


(28) 


for  n  =  m  and 

MA  A 
n  m 


M++  -  M+. 


M_  +  +  M_ 


ITT 

n  m 


+  o(6  6  ) 
n  nr 


for  n  i  m.  The  matrix  is  really  five  matrices.  One  is 

determined  using  values  in  the  above  equations,  and  then 
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the  others  by  through  in  turn.  The  error  terms  in  these 

equations  can  be  ignored  if  6  and  A  are  quite  small.  The 

n  m 

6  used  for  the  central  differences  technique  is 

A  =  DELTA  •  A  (29) 

n  n 

and  if  the  absolute  value  of  A^  is  larger  than  DELTA,  then 

6  =  DELTA  (30) 

n 

where  DELTA  is  another  small  positive  number.  In  this  prob¬ 
lem,  DELTA  was  initially  chosen  to  equal  1.0  x  10 

Selecting  Initial  Parameter 
Values 

Picking  nine  parameters  to  serve  as  an  initial  guess 

for  the  orbital  transfer  problem  is  far  from  a  simple  matter. 

Not  only  must  the  parameters  be  compatible,  but  they  should 

also  provide  an  initial  "miss"  relative  to  the  end  condition 

constraints  small  enough  to  allow  for  easy  convergence.  At 

this  point,  some  engineering  insight  was  applied  with  only 

marginal  success.  Since  the  actual  parking  orbit  has  an 

eccentricity  of  0.01497,  one  would  think  parameters  good  lot 

a  perfectly  circular  orbit  would  be  good  initial  guesses  for 

the  actual.  With  this  in  mind,  a  test  orbit  was  constructed 

in  an  F.CT  frame  with  ascending  node  on  the  X-axis  at  2.16737 
7 

x  10  feet,  inclined  22.8  degrees  at  a  velocity  of  25,618 
feet  per  second.  The  first  parameter,  T1 ,  was  determined  by 
calculating  the  time  of  flight  to  the  descending  node.  Since 
the  final  orbit  would  lie  in  the  X-Y  plane,  the  first  burn 
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was  aligned  to  eliminate  as  much  inclination  as  possible  using 


AV.  =  2  V  sin  9/2  (31) 

where  AV^  is  the  change  in  velocity  due  to  stage  burning,  V 
is  the  original  velocity  of  the  spacecraft  at  nodal  crossing, 
and  6  represents  the  angular  change  in  the  inclination.  This 
allowed  A1  and  B1  to  be  determined.  T2  was  chosen  to  be  the 
time  of  flight  to  the  next  nodal  crossing.  The  second  burn 
was  aligned  in  an  attempt  to  eliminate  the  remaining  inclina¬ 
tion  yielding  A2  and  B2.  T3  was  again  chosen  to  be  the  time 
to  the  next  nodal  crossing,  and  A3  and  B3  were  varied  to  find 
the  smallest  error  in  the  end  conditions.  This  orbit  and  the 
actual  orbit  differ  only  by  a  rotation  about  the  Z-axis  to 
align  the  nodes- -  ignoring  eccentricity.  This  angle  was  then 
added  to  the  three  alpha  parameters,  and  T1  was  adjusted  to 
reflect  the  coast  time  to  first  nodal  crossing  from  the 

initial  fix.  The  initial  miss  vector  that  this  procedure 

12 

yielded  was  of  magnitude  10  ,  largely  due  to  the  circular!  ty 

constraint.  Random  adjustment  of  individual  parameters  proved 
to  be  a  hopeless  means  of  reducing  the  error.  The  circular¬ 
ity  constraint  proved  to  be  far  too  sensitive  to  allow  for 
rapid  convergence.  To  decrease  the  sensitivity  of  this  con¬ 
straint,  the  circularity  was  multiplied  by  a  scaling  factor 
of  10  This  allowed  the  program  to  drive  down  the  errors 
in  the  other  four  constraints.  Each  time  this  occurred,  the 
scaling  factor  was  increased  by  a  factor  of  ten  and  the  pro¬ 
cess  repeated.  This  allowed  the  scaling  factor  to  eventually 
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be  removed  once  the  vicinity  of  the  functional  minimum  was 
reached. 
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IV.  Results 


Choosing  P,  Q,  and  DELTA 

The  selection  of  the  P  and  Q  scaling  factors  is  unique 
to  every  optimization  problem.  By  their  selection,  one  can 
select  several  options.  The  first  and  most  commonly  used  is 
to  choose  P  very  small  and  Q  large  as  compared  with  P.  By 
doing  this,  the  program  will  attempt  to  satisfy  the  end  con¬ 
straints  before  doing  any  optimization.  A  second  method  is 
to  choose  P  larger  than  Q  in  an  attempt  to  satisfy  the  end 
conditions  after  finding  an  optimal  region.  The  third  common 
method  is  to  choose  P  and  Q  equal  and,  thereby,  optimize 
while  trying  to  satisfy  the  end  constraints.  All  three 
methods  were  tried  on  this  problem,  with  only  the  first  being 
successful.  It  was  found  necessary  to  first  satisfy  the  end 
conditions  and  then  search  for  an  optimal  solution. 

Second-order  optimization  methods  are  very  unstable 
far  from  the  optimal  solution.  This  requires  very  small 
changes  in  the  control  parameters  through  each  iteration  to 
prevent  divergence.  The  magnitudes  of  these  "delta"  quanti¬ 
ties  can  be  controlled  using  P  and  Q.  Trial  and  error 
showed  that  for  P  equal  to  1.0  x  10  x  and  Q  equal  to  1.0 
x  10  6,  sufficiently  small  parameter  changes  were  generated 
to  allow  the  program  to  run  for  an  average  of  25  iterations 
before  divergence  occurred.  This  was  caused  by  the  numerical 


20 


routine  generating  an  excessively  large  change  in  one  of  the 
parameters.  The  program  was  modified  to  allow  continuous 
operation  by  examining  the  norm  of  these  parameter  changes. 
Initially,  whenever  the  norm  exceeded  0.38,  the  divergence 
was  seen  to  occur.  The  program  was  then  modified  to  reiniti¬ 
alize  all  the  Lagrange  multipliers  and  restart  itself  whenever 
this  was  observed.  As  the  program  began  to  approach  the 
vicinity  of  the  functional  minimum,  the  numerical  routine 
became  more  stable.  The  limit  on  the  norm  of  the  parameter 
changes  was  relaxed  to  0.5  and  finally  to  1.0.  This  simple 
procedure  proved  extremely  valuable  after  unsuccessfully 
attempting  to  use  a  straight  gradient  method. 

A  third  parameter,  DELTA,  also  showed  a  large  impact 
on  the  success  of  the  routine.  For  a  similar  problem  using 
differential  equations  of  motion,  DELTA  would  be  chosen  in 
the  10  *  to  10  ^  range.  This  insures  the  accurate  calcula¬ 
tion  of  the  and  matrices.  In  addition,  the  size  of 
the  parameter  changes  (da's)  is  directly  proportional  to 
DELTA,  since  the  "delta"  perturbing  values  equal  DELTA  times 
each  element  of  A.  In  this  problem,  DELTA  equal  to  10  6  was 
used  based  on  this  previous  experience.  It  was  found  that 
the  numerical  routine  stagnated  when  the  solution  was  being 
approached.  It  appears  the  function  contours  were  steep  to 
the  point  where  the  parameter  changes  equal  to  10  ^  times 
their  present  value  would  not  allow  the  program  to  move  deeper 
into  the  contour  valleys.  At  this  point,  several  iterations 
were  run  with  the  same  initial  parameter  values,  but  with 


different  DELTAS.  By  comparing  the  norms  of  the  constraint 
errors,  it  was  found  that  performance  improved  as  DELTA  was 
progressively  decreased.  Best  performance  resulted  from  a 
DELTA  of  10  which  was  then  used  to  find  the  solution  pre¬ 
sented  here. 

Gradient  Method 

The  gradient  method,  which  uses  first-order  information 
in  an  attempt  to  move  toward  an  optimal  solution,  proved  to 
be  a  dismal  failure.  Because  the  performance  index  was  de¬ 
fined  as  the  sum  of  three  times,  the  gradient  method  tried  to 
drive  all  three  to  zero.  Gradient  information  was  used, 
however,  to  calculate  initial  values  of  Lagrange  multipliers 
by  Eq  (18).  These  were  then  used  in  the  second-order  method. 

Numerical  Results 

The  second-order  method  used  in  the  computer  program, 
when  modified  to  run  continuously,  steadily  drove  toward  a 
solution  that  satisfied  the  end  constraints  and  then  toward 
an  optimal  solution.  The  parameters  that  satisfied  the  end 
constraints  changed  very  little  during  the  optimization  por¬ 
tion  of  the  program,  which  suggests  that  this  may  be  a  locally 
unique  solution.  The  optimization  iteration  process  was 
stopped  with  P  and  Q  equal  to  one  when  the  norm  of  the  para¬ 
meter  changes  was  less  than  10  ^  and  the  norm  of  the  Lagrange 
multiplier  changes  was  also  very  small.  The  A  matrix  solution 
and  the  sensitivity  matrix,  M^,  are  shown  in  the  following 


TABU-  HI 


Parameter  Solution  Set 


T1  =  2030.2449995  sec 

T2  = 

3847.461750268 

sec 

A1  =  -0.435361007  rad 

A2  = 

0.732341983313 

rad 

B1  =  0.7111815869  rad 

II 

rH 

PC 

2.731217478098 

rad 

T3  =  7470. 

83934000 

sec 

A3  =  1.30013074628 

rad 

B3  =  0.02236409011 

rad 

TABLE  TV 

End  Constraint  Errors 

M( 1)  =  zd  =  -2.5465  x  10‘10  ft/sec 

M( 2)  *  z  =  9.5367  x  10'7  ft 

M(3)  =  synchronous  speed  =  -1.0477  x  10 ft/sec 

M(4)  =  synchronous  distance  =  8.583  x  10  6  ft 

M(5)  =  circularity  (R  V)  =  1.318  x  10’2 

Total  Transfer  Time  =  13348.54608  sec  = 

3.70792  hrs 

Figures  3  through  8  depict  the  transfer  orbit  from 
several  points  of  view.  Figure  5  shows  a  result  that  is 
what  surprising.  The  transfer  orbit  remains  in  the  same 
as  the  parking  orbit.  In  other  words,  the  transfer  opts 
increase  the  altitude  and  velocity  oT  the  spacecraft  earl 
in  the  transfer  and  complete  the  inclination  change  in  th 
final  burn.  As  Fq  (31)  points  out,  the  cost  to  change 
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p  1  a  n  c 
t  o 

V 
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TABLE  V 


Position  and  Velocity  After  Firing  of  Each  Stage 


Position  and  Velocity  After 

Firing  First  Stage: 

x  =  -2.0372124  x  107  ft 

xd  =  1.11011967  x  104  ft/sec 

y  =  -7.6147903  x  106  ft 

yd  =  -2.21389677  x  104  ft/sec 

z  =  -1.3724214  x  106  ft 

zd  =  -9.3978524  x  103  ft/sec 

Position  and  Velocity  After 

Firing  Second  Stage: 

x  =  -2.885451299  x  106  ft 

xd  =  -2.05530541  x  104  ft/sec 

y  =  3.20133328  x  107  ft 

yd  =  -1.176502728  x  104ft/sec 

z  =  1.257707996  x  107  ft 

zd  =  -5.014671267  x  103ft/scc 

Final  Position  and  Velocity 

after  Third  Burn: 

x  =  1.379638  x  108  ft 

xd  =  464.30038  ft/sec 

y  =  -6.351478  x  106  ft 

yd  =  10085.31809  ft/sec 

z  -  -8.118152  x  1 0‘5  ft 

zd  =  2.8067  x  10"9  ft/sec 

inclination  is  less  as  the  spacecraft  velocity  decreases. 

The  point  at  which  geosynchronous  orbit  is  obtained  is  where 
the  spacecraft  velocity  is  lowest.  This  points  out  that  the 
assumption  used  to  ohtain  initial  values  for  the  parameters 
by  eliminating  the  inclination  of  the  transfer  orbit  was 
wrong.  Figure  6  depicts  the  entire  transfer  orbit  in  three 
dimensions.  Figure  7  is  the  same  as  Fig  6,  except  that  the 
scaling  on  the  Z-axis  has  been  changed  to  show  separation  of 
the  near  earth  portions.  Figure  8  is  the  confirmation  that 
the  final  orbit  is,  indeed,  circular  and  geosynchronous. 
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Fig  6.  Depiction  of  Transfer  Orbit  in  3-D 
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Fig  8.  Depiction  of  Final  Orbit  in  F.CI  X-Y  Plane 


The  matrix  which  is  the  partial  of  M  with  respect 
to  the  A  matrix  indicates  the  sensitivity  of  each  end  condi¬ 
tion  to  a  change  in  any  one  of  the  A  matrix  parameters. 


TABLP.  VT 


Matrix 


T1 

A1 

B1 

Ml 

9.9414 

-11045.0 

-17707.0 

M2 

-78377.0 

83614230.0 

147191047.0 

M3 

26.306 

-28542.0 

-45384.0 

M4 

-497674.0 

473346710.0 

756855918.0 

M5 

9.9097 

-7.73  x  1012 

-1.215  x  1013 

T2 

A2 

B2 

Ml 

-10.353 

3133.023 

1377.528 

M2 

-248.21 

4481256.0 

60289542.0 

M3 

-26.293 

7204.65 

6391.623 

M4 

304516.131 

-103134155.0 

-82096056.0 

M5 

1.596  x  109 

1.3178  X  1012 

3.077  x  1011 

T3 

A3 

B3 

Ml 

-0.802865 

0 

10995.048 

M2 

14068. 347 

0 

0 

M3 

-2.0329 

2449.40 

-  239.582 

M4 

50201.67 

0 

0 

M5 

-646330649.0 

-1.48  x  1012 

-7.567  x  109 

The  matrix  is  extremely  useful  when  one  is  interested 
in  how  changes  in  one  of  the  parameters  in  the  A  matrix  wi II 
effect  the  end  conditions.  As  an  example,  consider  the  last 
column,  which  represents  the  change  in  end  conditions  due  to 
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TABLE  VII 


Values  for  B3 


B3 

Ml 

10995.048 

M2 

0 

M3 

-239.582 

M4 

0 

M5 

-7.567  x  109 

a  change  in  B3.  These  numbers  indicate  that  if  the  third 
stage  was  misaligned  by  plus  one  radian,  it  would  increase 
the  zd  component  of  velocity  by  10995.048  feet  per  second. 
Since  B3  has  no  effect  on  position,  the  numbers  that  corres¬ 
pond  to  z  and  geosynchronous  distance  are  both  zero.  The 

satellite's  speed  would  decrease  by  240  feet  per  second, 

9 

and  the  circularity  of  the  orbit  by  7  x  1 0  .  This  clearly 
illustrates  the  use  of  the  matrix  as  well  as  the  wide 
differences  in  magnitudes  of  the  individual  errors. 
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V.  Application 

The  application  of  the  information  found  in  this  study 
is  both  simple  and  straightforward.  The  parking  orbit  and 
the  transfer  orbit  are  fixed  in  inertial  space  due  to  the 
choice  of  an  ECI  reference  frame.  The  earth,  on  the  other 
hand,  spins  relative  to  this  frame.  To  place  a  satellite  in 
geosynchronous  orbit  over  a  specific  spot  on  the  equator, 
one  has  basically  two  choices.  The  first  is  to  time  the 
launch  of  the  satellite  from  the  earth  into  its  parking  orbit 
and  then  along  the  transfer  orbit  such  that  the  final  posi¬ 
tion  of  the  transfer  orbit  coincides  with  the  desired  equa¬ 
torial  position.  A  second  and  more  practical  approach 
would  be  to  launch  into  the  parking  orbit  whenever  convenient. 
The  period  of  the  parking  orbit  is  roughly  90  minutes,  which 
allows  for  a  transfer  every  hour  and  a  half. 

Because  this  problem  was  worked  in  an  ECI  frame  of 
reference,  the  position  of  the  satellite  at  the  end  of  the 
transfer  relative  to  its  position  at  the  beginning  is  always 
constant.  Therefore,  the  angular  difference  in  their  posi¬ 
tions,  as  seen  in  the  X-Y  plane  and  measured  about  the  Z- 
axis,  is  constant.  By  allowing  for  the  rotation  of  the  earth, 
a  transfer  to  a  given  position  over  the  equator  can  be  calcul¬ 
ated  in  terms  of  longitude.  Optimal  launch  conditions  occur 
when  the  longitude  of  the  parking  orbit  initial  fix  is 
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114.90  degrees  greater  than  the  longitude  of  the  desired 
geosynchronous  orbital  position.  For  example,  if  the  initial 
orbit  was  obtained  at  longitude  23.90°E  and  the  transfer  begun, 
the  final  position  would  be  at  longitude  91.0°W,  which  corres¬ 
ponds  to  the  Galapagos  Island  chain,  after  a  flight  time  of 
3.70792  hours. 
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VI.  Recommendations 


One  of  the  major  assumptions  made  in  this  study  was 
that  the  three  engine  burns  were  impulsive.  As  was  shown, 
this  reduced  the  alpha  and  beta  controls  to  scalar  parameters. 
In  reality,  the  engines  operate  for  specific  lengths  of  time 
as  shown  in  Table  II.  Since  the  computer  program  has  the 
capability  to  determine  the  alpha  and  beta  controls  as  poly¬ 
nomial  functions  of  time,  it  would  be  interesting  to  remove 
the  impulsive  burn  assumption  and  rework  this  problem  using 
the  solution  found  here  as  a  good  initial  guess. 

This  study  found  only  one  solution  to  this  transfer 
problem  and  suggests  that  it  may  be  unique.  If  the  original 
parking  orbit  were  equatorial,  there  would  exist  two  mirror 
image  solutions.  It  would  be  interesting  to  see  if  the 
methods  used  in  this  study  would  yield  both  solutions. 

The  suboptimal  control  technique,  as  applied  to  this 
problem,  did  yield  a  valid  optimal  solution,  but  the  method 
did  prove  expensive  in  terms  of  computer  time.  It  would  be 
worthwhile  to  apply  one  or  two  other  second  order  or  quasi - 
second  order  optimization  methods  to  this  problem  to  see  if 
a  more  efficient  method  exists. 
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:l  00  -  PROGRAM  ORB  I T  (  INPUT  »  OUTPUT  ) 

110=  III  hr  NO  I  ON  DA<9>  >MAA3<9»9)  »MAA4<9r9)  >MAA5<9,9)  »WZ<9> 

120=  DIMENSION  A ( 9 ) » AA ( 9 ) »GA( 9 ) 

1  7,0  DIMENSION  M'S)  rMA(5»9>  »HAA1  (  9  *91  ,  MAA7  (9*9) 

1 4 0  *  0 1 M E N s 1 0 N  & < 6 >  * D B (6)  , B B  < A  ) , H P ( 6 ) , B PI <  6)  »  BP 2 (6 > i B P 3 (6)> H  P 4  <  A  > 

1  GO  =  Li  I  MENS  ION  BEL  <  9  )  ,  MAMAT  (5,5)  r  MAMAT 1(5*  5  )  *  MAC  AT  (  5  ) 

160=  DIMENSION  FAT ( 9 ) * FAA ( 9 , 9 ) * FAA I ( 9 * 9 > r MAFAAI ( 5 * 9 ) 

170  =  DIMENSION  W(5*5)  *WI(5*5)  *C(5)  *D(9)  ,  PNIJ  <  5  )  ,  DPNU  (  5  )  *PNUMA(9> 

180=  DIMENSION  E ( 9 ) *  FA ( 9 ) * FFAA ( 9 , 9 > * FFAAI < 9 * 9 ) 

190=  DIMENSION  GAA ( 9  *  9 ) 

200=  DIMENSION  MTXFXT ( 12 * 12 ) * FFAAB ( 22  *  22 ) 

210=  DIMENSION  MTX( 12* 12) *MTXT( 12*12) *CONA( 12) rCONB< 12) *MTXF( 12* 12) 

220=  COMMON/M  I SC/NPH *  NPI *  NA  * TF * DT *  N  *  S  * I.L » PP * TOT  (  8  )  *  FI  XT 

230=  COMMON  X  *  Y  *  2  *  XD»  YD *  ZD 

240=  COMMON/GUES/ ALPHA 

250=  REAL  M*  Ml  *  M2  *M3  *  M4  *  MM1 * MM2  *  MM3  *  MM 4  * MM5 * MM6 *  MM7 *  MM8  *  MA 

260=  REAL  M5  *  M6  *  M7  *  M8  *  M9  *  M10 

270=  REAL  MM9*MM10*MM11*MM12*MM13*MM14*MM15 

280=  REAL  MM16»  MM  17  *MM18*MM19*  MM20  *  MAA1  *MAA2*MAA3*MAA4*  MAA5 

290=  REAL  MAMAT  * MAMAT I *  MAGAT * MAFAA I *  NORMA *  NORMM r  NORMP 

300=C 

310=C  PROGRAM  SUBOPT  IS  A  SUBOPTIMAL  CONTROL  TECHNIQUE  USED  TO  FIND  AN 
320=0  APPROXIMATE  SOLUTION  TO  AN  OPTIMAL  CONTROL  PROBLEM.  THE  SOLUTION  TO 
330=C  MOST  OPTIMAL  CONTROL  PROBLEMS  IF  THE  CONTROLS  CANNOT  BE  SOLOED 
340=C  ANALYTICALLY  IS  TO  GUESS  THE  CONTROLS  AND  THE  LAGRANGE  MULTIPLIERS 
350=C  AND  SEE  IF  END  CONDITIONS  ARE  MET  AND  THE  PERFORMANCE  INDEX  MINIMIZED 
360=C  THE  SUBOPTIMAL  CONTROL  TECHNIQUE  ASSUMES  THE  CONTROLS  ARE  A  LINEAR 
370=C  COMBINATION  OF  POLYNOMIALS  WITH  UNKNOWN  COEFFICIENTS.  IN  SO  DOING 
380=C  THE  PROBLEM  IS  CHANGED  FROM  A  FUNCTIONAL  MINIMIZATION  PROBLEM  TO  A 
390=C  PARAMETER  OPTIMIZATION  PROBLEM. 

400=C  IN  THIS  EXAMPLE  THE  CONTROLS  ARE  NOT  FUNCTIONS  OF  TIME. 

4 1 0=C  THIS  SIMPLIFIES  THE  PROGRAM  TO  ONE  OF  PARAMETER  OPTIMIZATION. 

420=C 

430=C  THIS  VERSION  ADOPTS  THE  ORBITAL  TRANSFER  PROBLEM 
440=C 

450=C  A  -"MATRIX  OF  PARAMETERS ( T1 * A1 > B1 *  T2  *  A2  *  B2  *  T3  *  A3  *  B3 ) 

460--C  DA  ^MATRIX  OF  PARAMETER  CHANGES  (NP  X  1) 

470--C  NE  “NUMBER  OP  STATE  VARIABLES 

480=C  NP  “TOTAL  NUMBER  OF  PARAMETERS 

490= C  NC  “NUMBER  OF  END  CONDITIONS  TO  BE  SATISFIED 

500=C  P  “SCALING  FACTOR  PERFORMANCE  INDEX 

510=C  U  “SCALING  FACTOR  FOR  END  CONDITION  CONSTRAINTS  <M> 

520=C  G  =  SCALAR  PERFORMANCE  INDEX 

530 -C  GA  “  PARTIAL  DERIVATIVES  OF  G  WRT  A'S  (1  X  NP) 

540=C  GAA  *  SECOND  PARTIAL  DERIVATIVES  OF  G  WRT  A'S  (2  -  NP  X  NP) 

550=C  M  =  MATRIX  OF  PRESCRIBED  FINAL  CONDITIONS  ( NC  X  1  ) 

560=C  MA  =  PARTIAL  DERIVATIVES  OF  M  WRT  A'S  ( NC  X  NP) 

570=C  MAA  =  SECOND  PARTIAL  DERIVATIVES  OF  M  WRT  A'S  (2  -  NP  X  NP) 

500=C  B  “STATE  VARIABLES  ( X  *  Y *  Z  f  XD*  YD *. ZD  > 

590=C  E  =  AUGMENTED  PERFORMANCE  INDEX  <G  +  PNUT*M> 

600=C  FA  “  PARTIAL  DERIVATIVE  OF  F  WRT  A'S  <1  X  NP) 


610=C  FAT  =  PARTIAL  DERIVATIVES  OF  F  WRT  A'S  TRANSPOSED  (NP  X  1) 
620=C  FAA  =  SECOND  PARTIAL  DERIVATIVE  OF  F  WRT  A'S  <NP  X  NP) 
630=C  PNU  =  LAGRANGE  MULTIPLIERS  <NC  X  1) 

640=C  DPNl)=  CHANGE  IN  LAGRANGE  MULTIPLIERS  ( NC  X  1 ) 

650-C 

660=C  THE  CIRCULAR  PARKING  ORBIT  IS  INCLINED  22.8  DEGREES 

670=C  AT  AN  ALTITUDE  OF  2.1637E07  FT 

680=C 

690=  T  =0 . 0 

700-  X=1 .02931 2E07 

710=  Y= 1 ♦ 73235 4 E07 

720=  7=7.R81747E06 

'0=  XD=-2 . 248185E04 

,  -10=  YD=9 . 356206E03 

750=  ZD=7 .958385E03 

760=C 

770=C  THE  THREE  BURNS  ARE  ORIENTED  IN  TIME  AND  SPACE 

780=C  BY  NINE  PARAMETERS 

790=C 

800=  T 1=2030. 244999504 

0 1 0  =  A 1  -.4 3 5  3  6 100 7 0 0 5 5 

820=  Bl=. 71 11 815869911 

830=  12=3847.461750268 

840=  A2=. 7323419833131 

850=  B2=-2. 731217478098 

860=  13=7470.839340005 

870=  A3=l .300130746286 

880=  B3=. 0223640901 197 

890=C 

900=C  THE  FINAL  ORBIT  IS  CALCULATED  USING  THE  F  &  G  SERIES 
910=C  SOLUTION  TO  THE  TWO  BODY  PROBLEM. 

920=C 

930=  R=<X*X+Y*Y+Z*Z)**0.5 

940=  V= ( XD*XD+YD*YD+ZD*ZD ) **0 . 5 

950=  PRINT*. *X<0)=  '.X 

960=  PRINT*. 'Y<0>=  '.Y 

970=  PRINT*. *Z<0)=  '.Z 

980=  PRINT*. ■XD(0)=  *.XD 

990=  PRINT*, *YD(0>=  ".YD 

1000=  PRINT* » ' ZD  <  0 )  =  * . ZD 

1010=  PRINT*. 'RANGE*  '.R 

1020=  PRINT*. 'VELOCITY=  '.V 

1030=  PRINT*,'  ' 

1040=  NP=9 

1050=  NC=5 

1060=  NE=6 

1070=  ITER=0 

90=  MAX=100 

m <v90=  ZED=1 . 

1100=  DELTA*  1  .E--08 

1110=  IMET =2 

1120=  0=1. 

1130=  P=1 . 

1140=  FN0RMA=10.0 

1150=  FGN0RM=10.0 

1160-  CC 1  =  1. OE  -04 

1170=  CC2=1 .0E-02 

1100=  CC3= 1 • 0E-08 

1190=  CC4=1 . 0E-04 

1200-  CC5=1 . 0C-04 


1210=  B  <  1 ) =X 

1220=  B<2>=Y 

1230=  B  <  3 ) =Z 

1240=  B<4)=XD 

1250=  B ( 5 ) =YD 

1260=  B ( 6 ) =ZD 

1270=  A ( 1 ) =T1 

1280=  A ( 2 ) =A1 

1290=  A  <  3  )  =  El  1 

1300=  A<4 )=T2 

1310=  A  (  5 ) =A2 

20=  A  <  6 )  =E<2 

a ->30=  A  (  7  )  =T3 

1340=  A ( 8 ) =A3 

1350=  A  <  9 )  =F<3 

1360=555  ITER=0 

1370=  PNU(1)=~. 00001 

1380=  PNU (2 ) =8 . 0E-08 

1390=  PNU ( 3 ) =- . 0000537 

1400-  PNU  (  4  )  =  1 . 757 E  -07 

1410=  PNU ( 5 ) =-l . 7E~  1 1 

1420=  EiPNUd  )=DPNU(2)=0.0 

1430=  EiPNU  (  3  )  =  Ei PNU  ( 4 )  =DPNU  <  5 )  =0 . 0 

1 440  =  52  FORMAT ( 1 X , 5E 1 5 . 7 ) 

1450=50  FORMAT ( 1X,0E15.7) 

1460=C 

1470=C  DETERMINING  M  MATRIX  BY  INTEGRATING  DIFFERENTIAL  CONSTRAINTS 
1480=C 

1490=  DO  1  1=1, NP 

1500=1  DA<I)=0.0 

1510=1000  DO  2  1=1, NP 
1520=  A< I  )=A< I )  fDA( I ) 

1530=2  AA  < I ) =A  < I ) 

1540=  PRINT*  , *  * 

1550=  PRINT#* • ITERATION  NUMBER  *  *  ITER  *  *  P=  'rPr’  G=  ',0 

1560=  PRINT*,'  ' 

1570=  DO  3  1=1, NE 

1580=3  BB(I)=B(I> 

1590=  S=1.0 

1600=  CALL  TRANS <AA, BED 

1610=  M< 1 )=ZD 

1620=  M(2)=Z 

1630=  M<3)  =  (  <XEi*XD+YD*YEi+ZD*ZEi)**0.5)-l  .0096E04 

1640=  M  ( 4  >  =  < (X#X+Y*YfZ#Z)#*0.5)-l .3811E08 

1650=  M<5)=X*XB+Y*YD+Z*ZD 

1660=  M<5)=M(5)/ZED 

'  70=  G=AA ( 1 )  +AA  <  4 ) 1 AA ( 7 ) 

. ^80=  ITER=ITER+1 

1690=  PRINT*,*  • 

1700=  PRINT*, 'A  MATRIX' 

1710=  DO  666  1*1, NP 

1720=666  PRINT*, 'AC, I,')*  ',A(I) 

1730=  PRINT*,'  ' 

1740=  DO  4  1=1, NC 

l  750-  PNU (  I  )  =PNI I i  DTDPNIM  I  ) 

17o0  -A  PRINT*, *M( ' , I , ')=  *,M(I) 

1770=  NORMM- 0.0 

1780=  DO  5  1=1, NC 

1790  5  NORMM- NORMM 1M ( I )**2 

1800=  NORMM*SQRT (NORMM)  H Q 


1810=  PRINT* >•  * 

1820=  PRINT*> "NORM  OF  h  =  "fNORMMf*  P.I.  =  ’fG 

1830=C 

1840=C  DETERMINING  MA  AND  MAA  BY  CENTRAL  DIFFERENCES 
1850=C 

1860=  S=0.0 

1870=  DO  100  1=1 »NP 

1880=  DO  A  J=1 t NP 

1890=6  AA(J)=A(J) 

1900=  DEL ( I >=DELTA#A  < I ) 

1910=  IF <AE<3<DEL<  I )  >  .LE. DELTA)  DEL <  I  >  =DEL.TA 

20=  AA  < I ) =A  < I ) +DEL ( I ) 

1930=  DO  7  J=1 fNE 

1940=7  BP  <  J ) =B <  J) 

1950=  CALL  TRANS <AAf BP) 

1960=  M1=ZD 

1970=  M2=Z 

1980=  M3= ( <  XD*XD4 YD*YD+ZD*ZD>**0.5)-1 .0096E04 

1990=  M4=( (X*X» Y*YiZ*Z)**0.5>' 1 .3811E08 

2 000=  M 5 ■■■■■  X * X D i  Y * Y D  t  Z * Z D 

2010=  M5=M5/ZED 

2020=  G1=AA< 1 H AA<4> 4AA<7> 

2030=  AA  < I ) = A  < I ) -DEL ( I ) 

2040=  DO  8  J=  1 1 NE 

2050=0  BP ( J) =B ( J ) 

2060=  CALL  TRANS <AA» BP) 

2070=  M6=ZD 

2080=  M7=Z 

2090=  M8=< (XD*XD+YD*YD+ZD*ZD)**0.5>-1 .0096E04 

2100=  M9=  < (X*X+Y*Y+Z*Z)**0.5)-1.3811E08 

2110=  M10=X*XD+Y*YD»Z*ZD 

2120=  M10=M10/ZED 

2130=  G2=AA( 1 ) +AA<  4 ) +AA<  7 ) 

2140=  MA  < 1 f I )  =  (M1-M6 ) / <  2 . 0*DEL  < I )  ) 

2150=  MA  <  2  f I )  =  ( M2-M7 ) / <  2 . 0*DEL  < I  )  ) 

2160=  MA  <  3  f I )  =  < M3-M8  >/ <  2 . 0*DEL ( I ) ) 

2170=  MA  <  4 , 1 )  =  <  M4-M9 ) / <  2 . 0*DEL ( I ) ) 

2180=  MA  ( 5  f  I )  =  ( M5--M10 )  /  ( 2 . 0*DEL.  ( I ) ) 

2190=  GA< I )=<G1-G2)/(2,0#DEL< I ) ) 

2200=  DO  100  K=1fI 

2210=  IF  <  K . EO . I ) GO  TO  707 

2220=  IF< IMET.EO.l )G0  TO  100 

2230=  DO  9  J=1fNP 

2240=9  AA  <  J ) =A ( J ) 

2250=  AA  <  I  )  =A  <  I )  +DEL.  <  I ) 

2260=  DEL ( K ) =DELTA#A <  K ) 

70=  IF  ( AE<S  <  DEL  ( K )  )  ♦  LE .  DELTA )  DEL  <  K  )=DELTA 

2280=  AA <  K  >  =A < K) +DEL  <  K ) 

2290=  DO  10  J=1fNE 

2300=10  BP1(J)=B<J) 

2310=  CALL  TRANS (AAf BP 1 ) 

2320=  MM1=ZD 

2330=  MM2=Z 

23 10=  MM3=<  <XD*XD+YD*YD+ZD*ZD)**0.5)  --1  .0096E04 

2350  MM 4  =  f ( X*X+ Y*Y+7*Z ) *#0 . 5 ) - 1 .3811E08 

2360=  MM5=X*XP+Y*YD+Z*ZD 

2370=  MM5=MM5/ZED 

2380=  GG1 =AA< 1 ) +AA( 4 ) +AA<  7 ) 

2390=  DO  It  J=1fNP  4l 

2400=11  AA ( J) =A  <  J ) 


■ 


2410=  AA  < I >  =A  < I ) -f  DEL ( I ) 

2420=  AA ( K  )  =A  <  K ) -DEL ( K ) 

2430=  DO  12  J=1 , NE 

2440=12  BP2<  J>=B( J) 

2450=  CALL.  TRANS < AA , BP2 > 

2460=  MM6=ZD 

2470=  MM7=Z 

2480=  MM8=( <XD*XD+YIi*YD+ZD*ZD>**0.5)-l .0096E04 

2490=  MM9=<  <X*XTY*Y+Z*Z>**0.5)-1 .3811E08 

2500=  MM10=X*XD»Y*YD+Z*ZD 

2510=  MM 10=MM 10/ZED 

2520=  GG2=AA  < 1 ) -f  AA ( 4 ) +AA ( 7 ) 

’0=  DO  13  J=1 »  NP 

2^40=13  AA  <  J ) =A ( J ) 

2550=  AA ( I ) =A ( I ) -DEL ( I ) 

2560=  AA  <  K  >  =A  <  K ) +DEL ( K ) 

2570=  DO  14  J=1 > NE 

2580=14  BP3(J)=B(J) 

2590=  CALL  TRANS ( AA, BP3 ) 

2600=  MM11=7D 

2610  MM 12=7 

2620=  MM  13=  <  <XD#XD+YD#YD»ZD*ZD)**0.3)--1 .0096E04 

2630=  MM 14= ( ( X#X+Y*Y+Z*Z ) **0 . 5 ) - 1 . 381 1 E08 

2640=  MM15=X*XDiY*YD+Z#ZD 

2650=  MM15=MM15/7ED 

2660=  GG3=AA  ( 1. )  »  AA  <  4  )  +  AA  <  7  ) 

2670=  DO  15  J= 1 r  NP 

2680=15  AA ( J ) =A ( J ) 

2690=  AA ( I ) =A( I ) -DEL ( I  > 

2700=  AA( K ) =A<  K ) -DEL  <  K ) 

2710=  DO  16  J=1 ,  NE 

2720=16  BP4 ( J ) =B  <  J ) 

2730=  CALL  TRANS ( AA , BP4 > 

2740=  MM16=ZD 

2750=  MM17=Z 

2760=  MM18=<  <XD*XD+YD*YD+ZD*ZD)**0.5>-1 .0096E04 

2770=  MM19=< (X*X+Y*Y+Z*Z)**0.5>-1 .3811E08 

2780=  MM20=X*XD+Y*YD+Z*ZD 

2790=  MM20=MM20/ZED 

2800=  GG4=AA ( t >  f AA < 4 >  * AA < 7 > 

2810=  MAA1 < I,K)=MAA1 <K, I ) = < MM1 -MM6-MM 1 1 +MM16 >/ (4 .0*DEL < I >*DEL<K) ) 

2820=  MAA2  <  I » K ) =MAA2  <  K » I  >  =  <  MM2-MM7-MM1 2  f MM17 ) /<  4 . 0*DEL < I )#DEL ( K ) ) 

2830=  MAA3 <  I ,  K  ) =MAA3 ( K > I )  =  <  MM3-MM8-MM1 3+ MM18 ) / ( 4 . 0#DEL  < I > *DEL  <  K > ) 

2840=  MAA4  < I ,K > =MAA4 < K , I ) = < MM4-MM9-MM1 4 f MM1 9 ) / < 4 . 0*DEL < I >#DEL(K> ) 

2850=  MAA5  <  I  ,K>  =MAA5  <  K» I )  =  <  MM5-MM10-MM15+MM20  )  /  <  4  . 0*DEL  <  I )  #DEl.  ( K ) ) 

2860=  GAA  <  I ,  K  )=GAA  ( K » I )  =  ( GG1 -GG2 -GG3+ GG4  )  /  (  4 . 0*DEL.  ( 1  )  *DEL  <  N )  ) 

2870=  GO  TO  100 

00=707  MAAK I , I ) = < Ml -2 . 0*M < 1 )+M6)/<DEL( I >**2> 

«,,.j90=  MAA2  (  T  ,  I )  =  ( M2~2 . 0*M  ( 2 )  +M7  )  /  (  DEL  ( I  )  **2 ) 

2900=  MAA3(T . T ) = < M3-2 ♦ 0*M < 3 > +M8 ) / < DEL ( I >**2> 

2910=  MAA4( I , I )  =  < M4-2 . 0#M  <  4  >  TM9 ) / <  DEL < I >  *$2  > 

2920=  MAA5< I r I )=< M5-2 .0*M<5> EM 10 ) / < DEL < I ) **2 ) 

2930=  GAA  <I»I)=<Gl-2. 0*G+G2 ) / ( DEL  < I ) ##2 ) 

2940=100  CONTINUE 

2950=  IF  ( I MET- 1 ) 759 » 759 ,747 

2960  i 

2970-C  riNDING  INITIAL  LAGRANGE  MULTIPLIERS  AND  DA'S  (GRADIENT  rECII) 
2900  ■  C 

2990=759  DO  18  1=1, NC 

3000  DO  18  J=1,NC  l\;' 


•  f 


3010  =  MAMAT  < I »  J ) =0 . 0 

3020=  DO  10  K=1 »NP 

3030=18  MAMAT  <  I » J ) =MAMAT  <  I ,  J ) +MA < I , K ) *MA  <  J , K ) 

3040=  CALI  GAUSD  <  NC  »  1 . 0E-30 » MAMAT »  MAMAT  1 , 1  ER  »  NC  > 

3050=  DO  19  1=1, NC 

3060=  MAGAT ( I  )  =0 • 0 

3070=  DO  19  J=1,NP 

3080=19  MAGAT ( I ) =MAGAT ( I ) IMA <  I »  J  ) *GA ( J ) 

3090=  DO  20  I = 1 » NC 

3100=  PNU  < I ) =0 , 0 

3110=  DO  20  J=1  r  NC 

7  J  20  =  20  PNU  < I >  =PNU < I ) TMAMATI ( I ,  J ) * ( M ( J ) -MAGAT ( J  >  > 

...  30 --  PRINT*,  *  * 

3140=  DO  21  1=1, NC 

3150=21  PRINT*, *  PNU ( • ,1 , * )  =  ’ , PNU < I ) 

3160=  PRINT* , '  ‘ 

3170=  DO  34  1=1, NP 

3180=  E  < I ) =0 . 0 

3190=  DO  34  J=1 ,NC 

3200=3  !  Ed)  =E  <  I )  IMA  <  ,.J  *  I  )  *PNtJ  <  J ) 

3210-  DO  35  1=1, NP 

3220=  FA  < I >  =GA ( I )+E  < I ) 

3230=35  DA  < I ) =-P*FA ( I ) 

3240=  F  =  0 . 0 

3250=  DO  36  1=1 ,NP 

3260=36  F=F+DA< I )**2 

3270=  GRNORM=SQRT ( F ) 

3280=  DIFF= <  FGNORM-GRNORM ) /FGNORM 

3290=  FGNORM=GRNORM 

3300=  PRINT*,*  * 

3310=  PRINT*,* FA  MATRIX* 

3320=  PRINT  50 , < FA (I), 1  =  1, NP ) 

3330=  PRINT*,* 

3340=  PRINT*, 'GRADIENT  METHOD  DA'S* 

3350=  PRINT  50 , ( DA < I ) , I  =  1 , NP > 

3360=  PRINT*,*  * 

3370=  IFIDIFF.LT.CC2. AND. P.EQ. 1 ,0)  0=1.0 

3380=  IF ( V .EQ . 1 .0  >PRINT* , *  GRADIENT  METHOD  CONVERGENCE* 

3390=  PRINT*,*  * 

3400=  IF ( V. EG .1.0)  P= . 1 

3410=  ITER=0 

3420=  GO  TO  1000 

3430=C 
3440=C 
3450=C 
3460=C 
=A70=C 
80=C 
3490=C 
3500=C 
3510=C 

3520=C  CALCULATING  DPNU  AND  DA  (SECOND  ORDER  TECH) 

3530=C 

3540=747  DO  22  1=1, NP 

3550  PNUMA  < I  '  =0 . 0 

3560  DO  73  J=t,NC 

35 70  23  PNUMA ( I )=PNUMA (1)1 PNU ( J ) *MA ( J , I  ) 

3580=22  FAT  < I ) =GA < I ) +PNUMA ( I ) 

3590-  PRINT*.*  * 

3600=  PRINT*, *PNU'S*  4} 


36 10  = 
3620= 

3630  = 

3640  = 
3650= 
3660= 

3670  = 
3600=24 
3690  = 

3700  = 
3710=7007 
720  = 

3730  = 
3740= 

3750  = 

3760=25 

3770= 

3780  = 

3790  = 

3000  -■ 

3810=26 

3820= 

3830= 

3040= 

3350= 

3860=27 

3870= 

3880  = 
3890= 
3900=28 
3910= 

3920  = 

3930= 

3940= 

3950= 

3960= 

3970=29 

3980= 

3990  = 
4000= 
4010=30 
4020= 

4030  = 
4040=31 
4050= 
4060= 

170= 
4080= 
4090= 
4100= 
4110= 
4120= 
4130=32 
4  I  40  = 

4150 - 
4160  = 
4170= 

4180 = 

4190 

4200= 


PRINT  50 ,  ( PNU <I) , I  =  1 » NC ) 

PRINT*,  *  * 

PRINT*, 'FA  MATRIX* 

PRINT  50, (FAT(I) ,I=1,NP) 

DO  24  I = 1 , NP 
DO  24  J=1 , NP 

FAA < I , J ) =PNU ( 3 ) *MAA3 < I »  J ) +PNU <  4 ) *MAA4  < I , J ) +PNU ( 5 ) *MAA5 < I , J) 
FAA  < I , J  >=FAA  < I , J ) +PNU ( 1 ) *MAA1 ( I , J ) +PNU <  2 ) #MAA2  < I , J ) +GAA ( I , J ) 
PRINT*, *FAA=  * , FAA ( I , J ) 

NF =NP 

CALL  GAIJSD  (  NF  ,  1 . 0F-30  ,  FAA  ,  FAA  I  ,  JER  ,  9  ) 

DO  25  1=1, NC 
DO  25  J= 1 , NP 
MAFAAI < I , J ) =0 , 0 
DO  25  K= 1 , NP 

MAFAAI (I , J)=MAFAAI (I , J) fMA< I »K)*FAAI <K, J) 

DO  26  1=1, NC 
DO  26  J=1  ,  NC 
U< I, J) =0.0 
DO  26  K= 1 , NP 

W ( I ,  J ) =W ( I , J ) +MAFAA I ( I , K ) *MA  <  J , K ) 

CALL  GAIJSD  ( NC  ,  1 . 0E-30  ,  U ,  W 1 ,  KER ,  NC  ) 

DO  27  1=1, NC 
C < I > =0 . 0 
DO  27  J=1 ,NP 

C< I  )=C< I H MAFAAI ( I, J)*FAT( J) 

DO  28  1=1, NC 
riPNIK  I  >=0.0 
DO  28  J= 1 , NC 

DPNIJ  (  I  )  =DPNU  <I)+UI<I,J)*<  -P*C  (  J  )  +Q*M  (  J  )  ) 

PRINT*,*  * 

PRINT*, "DPNU'S* 

PRINT  50 , < DPNU (I),I=1,NC) 

DO  29  1=1, NP 
D  < I ) =0 . 0 
DO  29  J=1 , NC 

D ( I > -D ( I ) +MA <  J , I ) #DPNU ( J ) 

DO  30  1=1, NP 
DA< I) =0.0 
DO  30  J=1 » NP 

DA  < I ) =DA ( I ) +FAA I < I , J ) *  < -P#FAT (J)-D(J)) 

N0RMP=0.0 
DO  31  1=1, NC 
NORMP=NORNP+DPNU ( I ) **2 
NORMP=SBRT  <  NORMP ) 

PRINT*,*  * 

PRINT*, ’NORM  OF  DPNU'S=  *, NORMP 
PRINT*,*  * 

PRINT*,* DA  MATRIX* 

PRINT  50, <DA< I ) , 1=1 ,NP) 

N0RMA=0 ,0 
DO  32  1=1, NP 
NORMA=NORMA+DA ( I ) **2 
NORMA=SOPT ( NORMA ) 

DIF--  <  FNORMA-NORMA )  /FNORMA 
FNORMA=NORMA 
PRINT*,*  * 

PRINT*,* NORM  OF  DA'S  =  *, NORMA 
IF< ITER. GE« 10.0) GO  TO  555 
IF < NORMA. GE. 1 .0)00  TO  4371 
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4210=  0=1. 

4220=C 

4230=C  CONVERGENCE  CRITERIA 
4240=C 

4250=200  IFINORMM. LE.CC3. AND. NORMA. LE.CC4. AND. P.EQ. 1 .0)  GO  TO  201 

4260=  IF(P.GT.l.O)  P=1.0 

4270=  GO  TO  1000 

4200=4371  0=0/100. 

4290=  GO  TO  555 

4300=201  PRINT# >  *  * 

4310=  PRINT# »  *  ^CONVERGENCE#  PERFORMANCE  INDEX=  *  y  G  y  *  P=  • 

'20=757  STOP 

..>30=  END 

4340=C 
4350=C 

4360=C  MATRIX  INVERSION  SUBROUTINE 
4370=C 


4380=C 

4390= 

SUBROUTINE  GAUSD ( M  y EPS  y  D  y  C  »  KER *  LAY ) 

4400 

DIME:  NS  ION  B  <  LAY  y  LAY  )  y  C  <  LAY  y  LAY  )  ,  A  < 

4410  = 

DOUBLE  PRECISION  Z y A y X y S y RATIO y EP 

4420  = 

EP  =  EPS 

4430= 

N  =  M 

4440  = 

DO  100  J  =  1 y  N 

4450= 

DO  100  K  =  1 y  N 

4460= 

100 

A(JyK)  =  B  <  J  y  K ) 

4470= 

DO  1  1=1 yN 

4400  = 

DO  1  .1=1  yN 

4490  = 

1 

X (  I y  J )  =  O.ODO 

4500  = 

DO  2  K=1 yN 

4510= 

2 

X ( K , K )  =  l.ODO 

4520= 

10 

DO  34  L  =  1 y  N 

4530= 

KP=0 

4540= 

Z  =  O.ODO 

4550= 

DO  12  K=LyN 

4560= 

IF  <  Z-DABS  <  A  <  K  y  L ) ) )  Ilyl2yl2 

4570= 

11 

Z  =  DABS(A(KyL) ) 

4580= 

KP=K 

4590  = 

12 

CONTINUE 

4600= 

I F  (L.-KP  )13y20y20 

4610  = 

13 

DO  14  J=LyN 

4620= 

Z=A(Ly J) 

4630= 

A (L  y  J ) =A  <  KP  y  J ) 

4640= 

14 

A<KPy J)=Z 

4650= 

DO  15  J= 1 y  N 

4660= 

Z  =  X  <  L  y  J ) 

'70  = 

X  <  L  y  J ) =X  <  KP  y  J ) 

u80= 

15 

X  <  KP y  J ) =Z 

4690= 

20 

IF  < DABS  <  A ( L  y  L  >  > - EP  > 50 y  50  y 30 

4700= 

30 

IF (L-N)31 y34y34 

4710= 

31 

LP1=L+1 

4720= 

DO  36  K=LP1 y  N 

4730= 

IF(A(KyL>  )32y36y.32 

4740  = 

32 

RAT  10= A  <K  yl.  )  /A  (L  y  L.  > 

4  750 

DO  33  J  L P 1  y  N 

4760 

33 

A(  K y J ) =A ( K  y J ) - RATIO*A <  L  y  J ) 

4770 

DO  35  J= 1 y  N 

4700  = 

35 

X<Ky.l)=X(Ky  J)-RATIO*X<l..y  ,J) 

4790  = 

36 

CONTINItF 

4000= 

34 

CONTINUE 

4310" 

40 

no  43  1  =  1 ,  N 

4820= 

I I=N+1-I 

4830= 

DO  43  J=1 » N 

4840= 

s  =  o.ono 

4850= 

IF< II-N)41 ,43,43 

4860= 

41 

IIP1=II+1 

4870= 

no  42  K  =  I  IF'l  ,  N 

4880= 

42 

S=S+A< II»K>*X(K»J) 

4890  = 

43 

X(II»J)=(X(II,J)— S)/A(II,II)  ! 

4900  = 

KER=1  : 

4910= 

no  200  J  =  1,N 

4920= 

no  200  K  =  1 ,  N 

*930  = 

200 

C  ( J » K )  =  X  <  J ,  K ) 

740  = 

GO  TO  75  1 

4950= 

50 

KER=2  i 

4960= 

70 

F'RINT  71 

4970= 

71 

FORMAT ( 1X,*MATF<IX  SINGULAR  IN  GAUSD# ) 

4980= 

75 

CONTINUE 

4990= 

RETURN 

5000= 

END 

501 0=C 

5020=C 

5030=C 

5040  = 

SUBROUTINE  CHEBY  <  T ) 

5050= 

DIMENSION  7.(7) 

5060= 

COMMON/M ISC/NF'H ,  NP1 , NA  ,  TF  ,  N  r  S , l.L  , PP  ,  TOT  (  8  )  , FI  XT 

5070  = 

COMMON/QUES/ ALPHA 

5080= 

Z ( 1 ) =T 

5090= 

no  1  1=2,7 

5100= 

K  =  I~1 

5110=1 

Z< I >=Z(K>#T 

5120= 

TOT ( 1 ) =1 4  0 

5130  = 

TOT ( 2 ) =2 . 0#Z  < 1 ) - 1 . 0 

5140= 

T0T<3)=8.0*Z<2)-8.0*Z< 1 )+l ,0 

5150= 

T0T<4)=32.0*Z<3> -  48 . 0*Z  <  2  H18 . 0*Z  < 1 )-l .0 

5160= 

TOT  ( 5 )  =128 . 0*Z  <  4 )  -256 , 0*Z  <  3 )  +  160  ♦  0*Z  (  2 )  -32! .  0*Z  (1)41.0 

5170= 

TOT<6>=512.0*Z(5>-1280.0*Z(4)+1120.0*Z<3>-  400 . 0*Z ( 2 > +50 . 0*Z < 1 )-l .0 

5180= 

T0T<7)=2048.0*Z<6)-6144,0*Z(5>+6912.0*Z(4)~3584.0*Z(3H840.0*Z(2> 

5190= 

1  -72.0*Z(1)+1.0 

5200= 

T0T(8)=8192.0*Z<7>-28672.0#Z<6)+39424.0*Z(5)-26880.0*Z<4> 

5210= 

1  i9408.0*Z<3>-1568,0*Z<2>+98.0*Z< 1 )-l .0 

5220= 

RETURN  $  END 

5230= 

SUBROUT I NE  TRANS  <  A , B ) 

5240= 

DIMENSION  A < 9 ) , B < 6 ) 

5250= 

COMMON  X,Y»Z,XD»YD,ZB 

5260= 

X=B ( 1 ) 

5270  = 

Y=B (2 ) 

C:2Q0= 

Z=B (3) 

290= 

XD=B<4) 

5300= 

Yn=n<5) 

5310= 

ZD=B<6) 

5320  = 

T1=A<  1 ) 

5330= 

A1=A  <  2 ) 

5340= 

B1 =A  <  3 ) 

5350  = 

T2=A(4) 

5360= 

A2=A(5) 

5370 

B2=A  <  6 ) 

•.Kmjv 

T3=A  <  7  ) 

5390= 

A3=A<  8 ) 

5400  - 

B3=A ( 9 ) 

•  * 
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kii. . __ . 

_ 

_ _ _ _ _ J 

T0F=T1 

CALL  FNG  <  TOF » X » Y » Z » XD » YD » ZD ) 


5410= 

5420= 

5430=C 

5440=C  VELOCITY  CHANGE  DUE  TO  FIRING  OF  FIRST  STAGE 

5450=C 

5460=  IiXD=4242.175*C0S<Bl)*C0S<Al ) 

5470=  DYD=4242.175*CC]S<B1  >*SIN<A1  > 

5480=  DZD=4242«175*SIN<E<1> 

5490=  X»=XD+DXD 

5500=  YD=YD+DYD 

5510=  ZD=ZD+DZD 

5520=  T0F=T2 

30=  CALL  FNG ( TOF f X» Y » Z > XD » YD » ZD > 

j540=C 

5550=C  VELOCITY  CHANGE  DUE  TO  FIRING  OF  SECOND  STAGE 

5560=C 

5570=  nXD=9565.712*C0S(B2>*CQS<A2> 

5580=  DYD=9565 . 71 2*C0S ( B2 ) *S I N ( A2 ) 

5590=  DZD=9565 . 712*SIN ( B2 ) 

5600=  XD=XD<DXD 

56  :1.0  YD=YIH  DVD 

5620=  ZD=ZD+DZB 

5630=  T0F=T3 

5640=  CALL  FNG  <  TOF  *  X » Y  r Z  r XD » YD » ZD  > 

5650=C 

5660=C  VELOCITY  CHANGE  DUE  TO  FIRING  OF  THIRD  STAGE 

5670=C 

5680=  DXD= 1 0997 . 798*C0S ( D3 ) *COS  <  A3 ) 

5690=  DYD=10997.798*C0S(B3)*SIN<A3) 

5700=  DZD= 10997 • 798*SIN  < B3) 

5710=  XD-XD+DXD 

5720=  YD-YD+DYD 

5730=  ZD--ZD+DZD 

5740=  R=<X*X+Y*Y»Z*Z)**0.5 

5750=  V=  <  XD*XD+YD*YB4-ZD*ZD )  **0 . 5 

5760=  RETURN 

5770=  END 

5780=  SUBROUTINE  FNG ( TOF  » X  *  Y » Z  *  XD » YB  r  ZD  > 

5790=  AMU=1 . 4076468E16 

5800=  E=1 . 0E-06 

5810=  RQ=<X*X+Y*Y-»Z*Z)**0.5 

5820=  EPS® ( XD*XD+YD*YD+ZD*ZD ) /2 . O-AMU/RO 

5830=  A=-AMU/<2.0*EPS) 

5840=  XN=1 . 0 

5850=1  ZZ=XN*XN/A 

5860=  C=<1 .0-C0S<ZZ##0.5>  >/ZZ 

5870=  S=<  <ZZ**0.5>~SIN<ZZ**0,5) >/(ZZ#*l .5) 

100=  TN=(X*XD+Y*YD+Z*ZD)*C*<XN*XN)/<AMU**0.5> 

uB90~  TN=TN+< 1 .0-R0/A>*S*(XN**3.0> 

5900=  TN=TNfRO*XN 

5910=  TN=TN/<AMU*#0.5> 

5920=  RN=(XN*XN*C>+<X*XD+Y*YD+Z*ZD>*<1 .0-ZZ*S)/<AMU**0.5> 

5930=  RN=RN!RO# ( 1 . 0-ZZ*C ) 

5940=  DT  =  TOF-TN 

5950=  IF ( DT .LE .E ) GO  TO  99 

5960=  XN=XN+DT*  ( AMIJ**0 . 5  >  /RN 

5970=  GO  TO  1 

5900  99  F  =  1  .0  <XN*XN*r,)/<RO> 

5990=  0=T0F- <  XN**3 . 0 ) #5/ ( AMU**0 . 5 > 

6000  P=X 


•  4 
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6010=  Q=Y 

6020=  V=Z 

6030=  X=F*X»G*XD 

6040=  Y=F*Y+G*YD 

6050=  Z=F*ZfG*ZD 

6060=  GD=1 . 0  -  ( XN#XN )  *C/  (  ( X*X+Y*Y+Z*Z )  **0 . 5 ) 

6070=  FD=<AMU**0.5>*(ZZ*S-1  ,0>/(R0*<  ( X*X+Y*Y1Z*Z ) **0 . 5 )  ) 

<080=  FD=F0*XN 

090=  X0=F0*P+GD*XD 

6100=  YD«FD#G+GD* YD 

6110=  ZD=FD*U+GD*ZD 

6120=  RETURN 

6130=  END 

6140=*E0R 
6150=*E.0F 
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